Hormetic and synergistic effects of cancer treatments revealed by modelling combinations of radio - or chemotherapy with immunotherapy

Background Radio/chemotherapy and immune systems provide examples of hormesis, as tumours can be stimulated (or reduced) at low radio/chemical or antibody doses but inhibited (or stimulated) by high doses. Methods Interactions between effector cells, tumour cells and cytokines with pulsed radio/chemo-immunotherapy were modelled using a pulse differential system. Results Our results show that radio/chemotherapy (dose) response curves (RCRC) and/or immune response curves (IRC) or a combination of both, undergo homeostatic changes or catastrophic shifts revealing hormesis in many parameter regions. Some mixed response curves had multiple humps, posing challenges for interpretation of clinical trials and experimental design, due to a fuzzy region between an hormetic zone and the toxic threshold. Mixed response curves from two parameter bifurcation analyses demonstrated that low-dose radio/chemotherapy and strong immunotherapy counteract side-effects of radio/chemotherapy on effector cells and cytokines and stimulate effects of immunotherapy on tumour growth. The implications for clinical applications were confirmed by good fits to our model of RCRC and IRC data. Conclusions The combination of low-dose radio/chemotherapy and high-dose immunotherapy is very effective for many solid tumours. The net benefit and synergistic effect of combined therapy is conducive to the treatment and inhibition of tumour cells. Supplementary Information The online version contains supplementary material available at 10.1186/s12885-023-11542-6.


The model
Our model describing the interaction among the effector cells, tumor cells, and the cytokine (IL-2) is as follows 1−8 : T (t + ) = (1 − q 2 )T (t), where E(t) denotes the effector cells, which can activate immune-system cells such as cytotoxic T-cells, macrophages, and natural killer cells that are cytotoxic to the tumor cells.T (t) represents the tumor cells, and I L (t) denotes the concentration of cytokines including IL-2 in the single tumor-site compartment.
q 2 represents the tumor reduction rate in patients by exposure to radiotherapy and/or chemotherapy, and q 1 and q 3 depict side-effects of radiotherapy and/or chemotherapy on effector cells and the cytokine IL-2, respectively.Lastly, s 1 is a treatment term that represents an external source of effector cells such as LAK or TIL cells, and s 2 is a treatment term that represents an external input of cytokines into the system.The intrinsic growth rate of the tumor cells can be modelled by a linear growth term (in this case r 2 (T ) is constant) or as a type of limited-growth with a logistic or Gompertz type function.Here we chose the logistic growth function as follows: The dynamical evolution of the effector cell population is described by the first equation, where c represents the antigenicity of the tumor and µ 2 denotes the death rate of the effector cell population, the Michaelis-Menten term (i.e. the third term) reveals the saturated effects of an immune response stimulated by cytokines with a Michaelis-Menten constant g 1 and maximal response constant and I L (t + ) denote the values after the pulsed therapy at time t, where s 1 is an administration constant that represents an external source of effector cells each time; s 2 is an administration constant that represents an external input of IL-2 into the system each time, q 2 represents the tumor reduction rate in patients by combined exposure to radio/chemotherapy, and q 1 and q 3 depict side-effects of radio/chemotherapy on the effector cells and the IL-2 with 0 ≤ q 1 , q 2 < 1 and Note that the anti-tumor cytokines including IL-2 and IL-12 in cancer diseases have been measured in different cancer types, and the levels of IL-12 in sera were found to be increased in metastatic patients and decreased in gastrointestinal tumours, colon cancer and breast cancer after radio/chemotherapy 9,10 .
Those results indicate that radiotherapy and chemotherapy have dual effects on cytokines, i.e. in some cases they inhibit the growth of cytokines, and in other cases stimulate the growth of cytokines.Therefore, q 3 may be positive or negative and more generalized cases are considered in the main text to show the effects of random perturbations on the dynamics of model (S1.1).Moreover, 0 < q 3 < 1 indicates the reduction of levels of anti-tumor cytokines and q 3 < 0 means that the levels are increasing.

Existence and global stability of the tumor free periodic solution
The existence and stability of the tumor free periodic solution have been addressed in several publications 7,8 , but their authors only focused on the local stability of some special cases.So, in the following we present general results for the existence and stability of the tumor free periodic solution.In order to calculate the tumor free periodic solution for system (S1.1),we first let T = 0, and obtain the following subsystem Integrating the second equation dI L (t) dt = −µ 3 I L within the interval (nP, (n+ 1)P ], gives Thus, we have Let I * L be an equilibrium of the above equation, yielding and solving it with respect to I * L , we get Therefore, there exists a unique periodic solution for I L , denoted by I P L (t) with Further, integrating it within the interval (nP, (n + 1)P ], there is dt .
Then through some straightforward calculation we have Consequently, we have Note that, to ensure the existence of E * , we need As a conclusion, if condition (S2.4) holds true, then subsystem (S2.2) has a unique periodic solution with period P , denoted by (E P (t), I P L (t)).Here Therefore, we have the following proposition.
Proposition 1.If condition (S2.4) holds true, then subsystem (S2.2) has a unique periodic solution (E P (t), I P L (t)) with period P , which is globally stable.Proof.As the existence and uniqueness have been discussed above, here we only prove the global stability of the periodic solution (E P (t), I P L (t)).Firstly, we conduct the local stability analysis, which is determined by the two Floquet multipliers of the following matrix where Φ(P ) is a solution of the following equation at time point with Φ(0) = I (the identity matrix).Here, * indicates that the term is not necessary for the exact expression.
We can calculate the two Floquet multipliers of matrix M as follows: It is easy to find that λ 1 < 1 holds true in the condition (S2.4).This means that the periodic solution (E P (t), I P L (t)) is locally stable whenever it exists.As for the global stability, we can easily see that the component I P L (t) is globally stable and then by the theory of limitation systems we have the component E P (t) being also globally stable.This completes the proof.
Based on the above discussion, we know that model (S1.1) has a tumor free periodic solution (TFPS) (E P (t), 0, I P L (t)) in the condition (S2.4).Next, we determine the threshold condition for the stability of the tumor free periodic solution, which could be useful for designing a combination therapy strategy.
Similarly, the local stability of the TFPS is determined by the Floquet mul-tipliers of the following matrix and Φ(P ) is a solution of the following equation at time point with Φ(0) = I (the identity matrix).
Through some straightforward calculations, the three Floquet multipliers of matrix M are as follows: Similarly, we have that λ 1 < 1 holds true whenever the TFPS exists (i.e. in the condition (S2.4)).Thus, we can conclude that the TFPS is locally asymptoti- Next, we consider the global stability of the TFPS, that is, the global attraction of it.Consider the following system It is easy to verify that system (S2.5)has a periodic solution T P * (t) with period P , where for all t ∈ (nP, (n + 1)P ], which is globally asymptotical stable 11 .Further, there is dT /dt < r 2 (1 − bT )T .Thus, according to the comparison theorem of impulsive models, we have that there exists a t * and such that T (t) < T P * (t)+ when t > t * .Denote dt .
Assume λ * 2 < 1, then there exists an such that Definitely, λ * 2 < 1 implies λ 2 < 1.Furthermore, it is easy to see that dI L dt > −µ 3 I L .Thus, according to Theorem 1 and the comparison theorem of impulsive models, we have that there 1 .Thus, also according to Proposition 1 and the comparison theorem of impulsive models, we have that there exists a t * 2 such that Without loss of generality, we assume that the above inequality holds for all t ≥ 0. Hence, there is Therefore, we obtain = T (nP )δ.
Then, we verify that I L (t) → I P L (t) as t → ∞.Because T (t) → 0 as t → ∞, thus there exist and t * such that T (t) < and T (t) g3+T (t) < when t > t * .Here, we assume that T (t) < and T (t) g3+T (t) < hold true for all t > 0. Therefore, Then, by the comparison theorem of impulsive models, if µ 2 > p 1 , there exist an M and a t * such that E(t) < M when t > t * .
Without loss of generality, we also assume that E(t) < M holds for all t > 0 when µ 2 > p 1 .
Let Φ 1 = |I L (t) − I P L (t)| and taking the upper-right derivative of Φ 1 (t) yields It follows that Consequently, we have Φ 1 ((n + 1)P ) ≤ Φ 1 (nP + ) exp Thus, there is Φ 1 ((n + 1)P ) This means that Φ 1 ((n + 1)P ) → 0 as n → ∞ because can be small enough.That is, I L (t) → I P L (t) as t → ∞.Then, we prove that E(t) → E P (t) as t → ∞.Because I L → I P L and T (t) → 0 when t → ∞, there exists a t * such that I P L − < I L < I P L + and T (t) < for any positive number .Without loss of generality, we assume that Then, by considering where . As is small enough when t is big enough, we have Φ 2 ((n + 1)P ) → 0 as n → ∞, that is, E(t) → E P (t) as t → ∞.Based on the above discussion, we have the following conclusion.

Baseline parameter values and parameter sets used in the main text
In order to show the stability of the periodic solution (E P (t), I P L (t)) of subsystem (S2.2), we plot the solutions as shown in .Therefore, we focus on the parameter space at which the periodic solution (E P (t), I P L (t)) is globally stable, otherwise the E(t) will tend to infinity due to frequent administrations.The tumour-free periodic solutions and their stabilities for system (S1.1)have been shown in  The parameter set listed in Table S.1 of reference ( 1) is as follows: In order to increase the visibility of the figures shown in the main text, we make appropriate changes to get two parameter sets.For the first parameter set we let E 0 = T 0 = I L0 = 1 b , t s = r 2 and take The above confirms that the first parameter set used in the main text is as follows: x 0 = y 0 = z 0 = 1, c = 0.08, µ 2 = 0.1667, p 1 = 0.6917, g 1 = 10 For scaling all three variables we choose N 1 = 100, N 2 = 1000, N 3 = 1000, and we have the following second parameter set employed in the main text: c = 0.0128, µ 2 = 0.1667, p 1 = 0.6917, g 1 = 70 To fit the RCRC or IRC data sets mentioned in Section 2.4 of the main text, we got the two parameter sets listed in Table S.1.

The effects of checkpoints and treatment period
A realistic scenario is that the number of tumour cells is measured at the checkpoints when the patient's cancer was diagnosed and at subsequent followups, rather than being measured at the treatment point.To show this, without loss of generality we may assume that the patient has the disease checked periodically with a period P 1 (P 1 denotes the monitoring period), and the combination therapy takes place within two checkpoints nP 1 and (n + 1) This means that there exists a positive real number P satisfying 0 ≤ P ≤ P 1 , and at each time point nP 1 + P the combination therapy is applied.Therefore, P ], we get Hence, we have I L ((nP 1 +P ) + ) = (1−q 3 )I(nP 1 +P )+s 2 = (1−q 3 )I(nP 1 )e −µ3P + s 2 .Further, integrating the equation dI L (t) dt = −µ 3 I L in the interval (nP 1 + P, (n + 1)P 1 ] with an initial condition I((nP 1 + P ) + ), we obtain I((n+1)P 1 ) = I((nP 1 +P ) + )e −µ3(P1−P ) = ((1−q 3 )I(nP 1 )e −µ3P +s 2 )e −µ3(P1−P ) .
Let I((n + 1)P 1 ) = I(nP 1 ) and solving the above equation, we can get a fixed point I L * , which is given by Because there is no implementation of the control strategy at the time point nP 1 , we can also write the periodic solution of I L (t) in the following form: nP1) , t ∈ (nP 1 + P, (n + 1)P 1 + P ]. (S4.9) Note that, here I * * L is equal to the I * L of the former model (i.e.model (S2.2)).Therefore the periodic solutions of the two models in terms of I L are the same if we take the pulse periods as the same.Further, also because there is no control at nP 1 , we can consider the existence of the periodic solution of E(t) in the interval (nP 1 +P, (n+1)P 1 +P ].As the periodic solution of I L is the same, we can obtain the same periodic solution of E(t) by substituting the periodic solution of I L into model (S2.2).All these results confirm that if we fix the treatment period, then the treatment time point does not affect the existence and stability of the tumour free periodic solution.Thus, the interesting question is how does the number of tumour cells measured at checkpoints during the monitoring of a patient after diagnosis rather than at the treatment time influence the dynamics of tumour cells?To address this, we carry out some numerical investigations about this in the following.huge difference in monitoring or detecting tumour cells at different times will have a great impact on the treatment plan, and then on the treatment effect.Therefore, it is very important to design a reasonable treatment and detection scheme for the treatment of tumours.5.5556, µ 3 = 55.556,p 2 = 27.7778,g 3 = 1, q 1 = 0.08, s 1 = 0.5, q 2 = 0.2, q 3 = 0.08, s 2 = 0.5, P 1 = 12 in (C) and (D), P 1 = 17 in (E) and (F).

p 1 .
The growth rate of the tumor cell population is given by the second equation with an intrinsic growth rate r 2 and a carrying capacity parameter b, the second term also shows the saturated effects of an immune response on the tumor cells stimulated by effector cells with two constants g 2 and a.The changing rate of the cytokine concentration is given by the third equation, which shows the effects of interactions between effector cells and tumor cells on the growth of cytokines through Michaelis-Menten kinetics with two constants g 3 and p 2 , and µ 3 denotes the degradation rate of the cytokines.See reference(1) for explanations of all of the parameters in more detail.The second part, consisting of three impulsive maps, describes the effects of pulses of radio/chemo-immunotherapy applied at period P (n = 1, 2, • • • ) on the effector cells, tumor cells and IL-2.E(t + ), T (t + )
If we fix the parameter values as those shown in Fig.S.3 and let P increase from 3 to 17, then the time series of E(t) and T (t) have been generated in each subplot.Although the whole solution curves seem to be similar, the phase positions of the extreme points of each solution curve change greatly as the time points P of combined treatment vary.To reveal the effects in more detail, we plot various outputs in Fig.S.4 for the different combined treatment points, i.e.P = 5 or 11.The time series have been shown in Fig.S.4(A) and (B); The number of effector cells and tumour cells at the each checkpoint nP 1 are shown in Fig.S.4(C), and the number of effector cells and tumour cells at the each monitoring point nP 1 + P are shown in Fig.S.4(D).While the relations between E n and E n+1 , T n and T n+1 , E P n and E P n+1 , T P n and T P n+1 are shown in Fig.S.4(E) and (F) for P = 5, from which we can see that the differences are obvious including the maximum amplitudes and iteration patterns.Similar outputs are shown in Fig.S.4(G-J) for P = 11.The difference in the number of tumour cells between the monitoring point and the treatment point can also be discussed through bifurcation analyses, as shown in Fig.S.5 for one parameter bifurcation diagrams with respect to q 2 , and in Fig.S.6 for two-parameter bifurcation diagrams with respect to q 2 and P .One parameter bifurcation diagrams reveal that the numbers of tumour cells in the monitoring site and the treatment site are very different, and the bi-stable regions also change slightly, as shown in Fig.S.5(K) and (L), which indicate that the RCRCs can be significantly affected if we record the number of tumour cells at different time points.But we can see that the values of T n (i.e.measures of the number of tumour cells before the treatment) do not change too much as the treatment period, P , changes.However, the two-parameter bifurcation diagrams shown in Fig.S.6 reveal that the numbers of tumour cells detected at two different time points may change significantly with the change of parameters.There is no doubt that the

Table S .
1: The parameter sets for fitting the RCRC or IRC data sets.